/* recmul.cpp: bcmath library file. */
/*
    Copyright (C) 1991, 1992, 1993, 1994, 1997 Free Software Foundation, Inc.
    Copyright (C) 2000 Philip A. Nelson
    Copyright (C) 2013 XChinux<XChinux@163.com>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.  (COPYING.LIB)

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to:

      The Free Software Foundation, Inc.
      59 Temple Place, Suite 330
      Boston, MA 02111-1307 USA.

    You may contact the author by:
       e-mail:  philnelson@acm.org
      us-mail:  Philip A. Nelson
                Computer Science Department, 9062
                Western Washington University
                Bellingham, WA 98226-9062
       
*************************************************************************/

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdarg.h>
#include "bcmath_p.h"

/* Recursive vs non-recursive multiply crossover ranges. */
#if defined(MULDIGITS)
#include "muldigits.h"
#else
#define MUL_BASE_DIGITS 80
#endif

int mul_base_digits = MUL_BASE_DIGITS;
#define MUL_SMALL_DIGITS mul_base_digits/4

/* Multiply utility routines */

static bc_num new_sub_num(int length, int scale, char *value)
{
    bc_num temp;
    
    temp = (bc_num)malloc(sizeof(bc_struct));
    temp->n_sign = PLUS;
    temp->n_len = length;
    temp->n_scale = scale;
    temp->n_refs = 1;
    temp->n_ptr = NULL;
    temp->n_value = value;
    return temp;
}

static void _bc_simp_mul(bc_num n1, int n1len, bc_num n2, int n2len, 
        bc_num *prod, int /*full_scale*/)
{
    char *n1ptr, *n2ptr, *pvptr;
    char *n1end, *n2end;		/* To the end of n1 and n2. */
    int indx, sum, prodlen;

    prodlen = n1len+n2len + 1;

    *prod = bc_struct::new_num(prodlen, 0);

    n1end = (char *) (n1->n_value + n1len - 1);
    n2end = (char *) (n2->n_value + n2len - 1);
    pvptr = (char *) ((*prod)->n_value + prodlen - 1);
    sum = 0;

    /* Here is the loop... */
    for (indx = 0; indx < prodlen - 1; indx++)
    {
        n1ptr = (char *) (n1end - MAX(0, indx-n2len+1));
        n2ptr = (char *) (n2end - MIN(indx, n2len-1));
        while ((n1ptr >= n1->n_value) && (n2ptr <= n2end))
            sum += *n1ptr-- * *n2ptr++;
        *pvptr-- = sum % BASE;
        sum = sum / BASE;
    }
    *pvptr = sum;
}


/* A special adder/subtractor for the recursive divide and conquer
   multiply algorithm.  Note: if sub is called, accum must
   be larger that what is being subtracted.  Also, accum and val
   must have n_scale = 0.  (e.g. they must look like integers. *) */
static void _bc_shift_addsub(bc_num accum, bc_num val, int shift, int sub)
{
    signed char *accp, *valp;
    int  count, carry;

    count = val->n_len;
    if (val->n_value[0] == 0)
        count--;
    assert(accum->n_len+accum->n_scale >= shift + count);
  
    /* Set up pointers and others */
    accp = (signed char *)(accum->n_value +
            accum->n_len + accum->n_scale - shift - 1);
    valp = (signed char *)(val->n_value + val->n_len - 1);
    carry = 0;

    if (sub)
    {
        /* Subtraction, carry is really borrow. */
        while (count--)
        {
            *accp -= *valp-- + carry;
            if (*accp < 0)
            {
                carry = 1;
                *accp-- += BASE;
            }
            else
            {
                carry = 0;
                accp--;
            }
        }
        while (carry)
        {
            *accp -= carry;
            if (*accp < 0)
                *accp-- += BASE;
            else
                carry = 0;
        }
    }
    else
    {
        /* Addition */
        while (count--)
        {
            *accp += *valp-- + carry;
            if (*accp > (BASE-1))
            {
                carry = 1;
                *accp-- -= BASE;
            }
            else
            {
                carry = 0;
                accp--;
            }
        }
        while (carry)
        {
            *accp += carry;
            if (*accp > (BASE-1))
                *accp-- -= BASE;
            else
                carry = 0;
        }
    }
}

/* Recursive divide and conquer multiply algorithm.  
   Based on 
   Let u = u0 + u1*(b^n)
   Let v = v0 + v1*(b^n)
   Then uv = (B^2n+B^n)*u1*v1 + B^n*(u1-u0)*(v0-v1) + (B^n+1)*u0*v0

   B is the base of storage, number of digits in u1,u0 close to equal.
*/
static void _bc_rec_mul(bc_num u, int ulen, bc_num v, int vlen, bc_num *prod,
	     int full_scale)
{ 
    bc_num u0, u1, v0, v1;
    //int u0len, v0len;
    bc_num m1, m2, m3, d1, d2;
    int n, prodlen, m1zero;
    int d1len, d2len;

    /* Base case? */
    if ((ulen+vlen) < mul_base_digits
      || ulen < MUL_SMALL_DIGITS
      || vlen < MUL_SMALL_DIGITS )
    {
        _bc_simp_mul (u, ulen, v, vlen, prod, full_scale);
        return;
    }

    /* Calculate n -- the u and v split point in digits. */
    n = (MAX(ulen, vlen)+1) / 2;

    /* Split u and v. */
    if (ulen < n)
    {
        u1 = bc_struct::copy_num(BCG(_zero_));
        u0 = new_sub_num (ulen,0, u->n_value);
    }
    else
    {
        u1 = new_sub_num (ulen-n, 0, u->n_value);
        u0 = new_sub_num (n, 0, u->n_value+ulen-n);
    }
    if (vlen < n)
    {
        v1 = bc_struct::copy_num(BCG(_zero_));
        v0 = new_sub_num (vlen,0, v->n_value);
    }
    else
    {
        v1 = new_sub_num (vlen-n, 0, v->n_value);
        v0 = new_sub_num (n, 0, v->n_value+vlen-n);
    }
    bc_struct::_bc_rm_leading_zeros (u1);
    bc_struct::_bc_rm_leading_zeros (u0);
    //u0len = u0->n_len;
    bc_struct::_bc_rm_leading_zeros (v1);
    bc_struct::_bc_rm_leading_zeros (v0);
    //v0len = v0->n_len;

    m1zero = u1->is_zero() || v1->is_zero();

    /* Calculate sub results ... */

    bc_struct::init_num(&d1);
    bc_struct::init_num(&d2);
    bc_struct::sub(u1, u0, &d1, 0);
    d1len = d1->n_len;
    bc_struct::sub(v0, v1, &d2, 0);
    d2len = d2->n_len;


    /* Do recursive multiplies and shifted adds. */
    if (m1zero)
        m1 = bc_struct::copy_num(BCG(_zero_));
    else
        _bc_rec_mul (u1, u1->n_len, v1, v1->n_len, &m1, 0);

    if (d1->is_zero() || d2->is_zero())
        m2 = bc_struct::copy_num(BCG(_zero_));
    else
        _bc_rec_mul (d1, d1len, d2, d2len, &m2, 0);

    if (u0->is_zero() || v0->is_zero())
        m3 = bc_struct::copy_num(BCG(_zero_));
    else
        _bc_rec_mul (u0, u0->n_len, v0, v0->n_len, &m3, 0);

    /* Initialize product */
    prodlen = ulen + vlen + 1;
    *prod = bc_struct::new_num(prodlen, 0);

    if (!m1zero)
    {
        _bc_shift_addsub(*prod, m1, 2 * n, 0);
        _bc_shift_addsub(*prod, m1, n, 0);
    }
    _bc_shift_addsub (*prod, m3, n, 0);
    _bc_shift_addsub (*prod, m3, 0, 0);
    _bc_shift_addsub (*prod, m2, n, d1->n_sign != d2->n_sign);

    /* Now clean up! */
    bc_struct::free_num(&u1);
    bc_struct::free_num(&u0);
    bc_struct::free_num(&v1);
    bc_struct::free_num(&m1);
    bc_struct::free_num(&v0);
    bc_struct::free_num(&m2);
    bc_struct::free_num(&m3);
    bc_struct::free_num(&d1);
    bc_struct::free_num(&d2);
}

/* The multiply routine.  N2 times N1 is put int PROD with the scale of
   the result being MIN(N2 scale+N1 scale, MAX (SCALE, N2 scale, N1 scale)).
   */
void bc_struct::multiply(bc_num n1, bc_num n2, bc_num *prod, int scale)
{
    bc_num pval; 
    int len1, len2;
    int full_scale, prod_scale;

    /* Initialize things. */
    len1 = n1->n_len + n1->n_scale;
    len2 = n2->n_len + n2->n_scale;
    full_scale = n1->n_scale + n2->n_scale;
    prod_scale = MIN(full_scale,MAX(scale,MAX(n1->n_scale,n2->n_scale)));

    /* Do the multiply */
    _bc_rec_mul (n1, len1, n2, len2, &pval, full_scale);

    /* Assign to prod and clean up the number. */
    pval->n_sign = ( n1->n_sign == n2->n_sign ? PLUS : MINUS );
    pval->n_value = pval->n_ptr;
    pval->n_len = len2 + len1 + 1 - full_scale;
    pval->n_scale = prod_scale;
    _bc_rm_leading_zeros (pval);
    if (pval->is_zero())
        pval->n_sign = PLUS;
    bc_struct::free_num(prod);
    *prod = pval;
}
